##########################################################################################

library(data.table)
library(optparse)
library(parallel)
library(dplyr)
suppressWarnings(library(edgeR, quietly = T))

##########################################################################################

option_list <- list(
    make_option(c("--sample_list_file"), type = "character"),
    make_option(c("--rsem_file"), type = "character"),
    make_option(c("--gtf_file"), type = "character"),
    make_option(c("--out_path"), type = "character")
)

if(1!=1){
    
    sample_list_file <- "~/20220915_gastric_multiple/rna_combine/analysis/config/tumor_normal.list"
    rsem_file <- "~/20220915_gastric_multiple/rna_combine/analysis/RSEM/CombineCounts.FilterLowExpression.tsv"
    out_path <- "~/20220915_gastric_multiple/rna_combine/analysis/images/TMM"
    gtf_file <- "~/ref/GTF/gencode.v19.ensg_genename.txt"

}

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

sample_list_file <- opt$sample_list_file
rsem_file <- opt$rsem_file
out_path <- opt$out_path
gtf_file <- opt$gtf_file

dir.create(out_path , recursive = T)

##########################################################################################

info <- data.frame(fread(sample_list_file))
dat_tpm <- data.frame(fread(rsem_file))
dat_gtf <- data.frame(fread(gtf_file , header = F))
colnames(dat_gtf) <- c("gene_id" , "Hugo_Symbol")

##########################################################################################

class_type <- c( "Normal" , "IM" , "IGC" , "DGC")

##########################################################################################

dat_tpm_all <- dat_tpm[,-2]
colnames(dat_tpm_all) <- gsub( "[.]" , "-" ,  colnames(dat_tpm_all) )

##########################################################################################
## 构建分组信息
condition <- sapply(strsplit(colnames(dat_tpm_all)[-1] , "_") , "[" , 2)
condition <- factor(sapply(strsplit(condition , "-") , "[" , 1))    

#individual <- factor(sapply(strsplit(colnames(dat_tpm_all)[-1] , "_") , "[" , 1))
#coldata <- data.frame(row.names = colnames(dat_tpm_all)[-1] , condition , individual)

## counts变为整数
readCount <- round(dat_tpm_all[,-1])
rownames(readCount) <- dat_tpm_all$gene_id

##########################################################################################
## https://rpubs.com/LiYumei/806213
conditions <- condition

## Count matrix preprocessing using edgeR package
y <- DGEList(counts = readCount , group = conditions)

## 去除TPM低表达的基因
## Remove rows consistently have zero or very low counts
#keep <- filterByExpr(y)
#y <- y[keep,keep.lib.sizes=FALSE]

## Perform TMM normalization and transfer to CPM (Counts Per Million)
## 统一标化
y <- calcNormFactors(y,method="TMM")
count_norm <- cpm(y)
count_norm <- as.data.frame(count_norm)

##########################################################################################
## 标准化的矩阵,用于后续样本间的表达比较
## 后续除软件要求计算评分使用TPM，否则均使用该矩阵
normalized_counts <- count_norm
normalized_counts <- data.frame(normalized_counts)
normalized_counts$gene_id <- rownames(normalized_counts)
colnames(normalized_counts) <- gsub( "[.]" , "-" ,  colnames(normalized_counts) )
image_name <- paste0( out_path , "/CombineCounts.FilterLowExpression.TMM.tsv" )
write.table( normalized_counts , image_name , row.names = F , sep = "\t" , quote = F )